function [EQD2,Effect_eq,n_0,n_f]=EQD2_calculator(Target_effect,par_long,dep_a,dep_b)
%This function calculates EQD2, total dose of a issoeffective treatment schedule with
%2 Gy/fraction using a SF expression with functional dependencies dep_a and
%dep_b. A bisection method is employed to find the number of 2 Gy
%radiotherapy fractions required to achieve the target effect.

%Input: 
%    1. Target effect: the isoeffect to achieve with a 2Gy/fraction fractionation schedule
%    2. par_long: TCP input parameter values
%    3. dep_a: functional dependence with dose of the alpha term in the expression of the surviving fraction
%    4. dep_b: functional dependence with dose of the alpha term in the expression of the surviving fraction


%Output: EQD2: Isoeffective total dose in 2Gy/fraction fractionation schedule


%TCP model parameter values:
alphabeta=par_long(1);
a=par_long(2);
b=par_long(3);
lambda=par_long(6);
tk=par_long(7);
alpha=par_long(8);
%Value= Effect corresponding to n fractions of 2 Gy
n0=1:200;
%EQD2 calculation for LQL model (lines 21 to 60)
    if dep_a~=-1
        vector_2Gy_effects=2*n0*alpha.*(1+a.*2.^dep_a+2.*(1+b.*2.^dep_b)/alphabeta)-lambda*max(0,7*floor(n0/5)+(n0/5-floor(n0/5))*5-2*double((n0/5-floor(n0/5))==0)-(n0 - fix(n0))-double(mod(n0,1) == 0)-tk);
        %Calculate EQD2 if all cohorts have an issoeffective treatment with 2Gy/#
        if (find(vector_2Gy_effects>max(Target_effect),1,'first'))>0    
            for i=1:length(Target_effect)
                n=find(vector_2Gy_effects>max(Target_effect(i)),1,'first');
                n_0(i)=n;
                value=2*n*alpha.*(1+a.*2.^dep_a+2.*(1+b.*2.^dep_b)/alphabeta)-lambda*max(0,7*floor(n/5)+(n/5-floor(n/5))*5-2*double((n/5-floor(n/5))==0)-(n - fix(n))-double(mod(n,1) == 0)-tk);
                delta=1;
                while (value <0.9999*Target_effect(i) || value>1.0001*Target_effect(i))
                    if value > 1.0001*Target_effect(i)
                        n=n-delta;
                        t_trat=7*floor(n/5)+(n/5-floor(n/5))*5-2*double((n/5-floor(n/5))==0)-(n - fix(n))-double(mod(n,1) == 0); 
                        value_integern=floor(n)*2*alpha.*(1+a.*2.^dep_a+2.*(1+b.*2.^dep_b)/alphabeta);
                        d2=2*(n - fix(n));%take the fractional of n, multiply by 2 to get the last fraction dose
                        value_fractionaln=d2*alpha.*(1+a.*d2.^dep_a+d2.*(1+b.*d2.^dep_b)/alphabeta);
                        value=value_integern+value_fractionaln-lambda*max(0,t_trat-tk);
                        delta=delta/2;
                    elseif value < 0.9999*Target_effect(i)
                        n=n+delta;
                        t_trat=7*floor(n/5)+(n/5-floor(n/5))*5-2*double((n/5-floor(n/5))==0)-(n - fix(n))-double(mod(n,1) == 0); 
                        value_integern=floor(n)*2*alpha.*(1+a.*2.^dep_a+2.*(1+b.*2.^dep_b)/alphabeta);
                        d2=2*(n - fix(n));
                        value_fractionaln=d2*alpha.*(1+a.*d2.^dep_a+d2.*(1+b.*d2.^dep_b)/alphabeta);
                        value=value_integern+value_fractionaln-lambda*max(0,t_trat-tk);
                        delta=delta/2;
                    end          
                end
                Effect_eq(i)=value;
                n_f(i)=n;
                EQD2(i)=2*n;
            end                        
        else
    %If there is no issoeffective treatment with 2 Gy/# return NaN
            EQD2=NaN;
            Effect_eq=NaN;
            n_0=NaN;
            n_f=NaN;               
        end    
    else  
%EQD2 calculation for LQ or modified LQ models (lines 63 to 120)      
        vector_2Gy_effects=n0.*alpha.*(2+2.*(a.*2+exp(-a.*2)-1)./(alphabeta.*a^2))-lambda*max(0,7*floor(n0/5)+(n0/5-floor(n0/5))*5-2*double((n0/5-floor(n0/5))==0)-(n0 - fix(n0))-double(mod(n0,1) == 0)-tk);
        %Calculate EQD2 if all cohorts have an issoeffective treatment with 2Gy/#
        if (find(vector_2Gy_effects>max(Target_effect),1,'first'))>0
            for i=1:length(Target_effect)
                n=find(vector_2Gy_effects>max(Target_effect(i)),1,'first');
                n_0(i)=n;
                value=n.*alpha.*(2+2.*(a.*2+exp(-a.*2)-1)./(alphabeta.*a^2))-lambda*max(0,7*floor(n/5)+(n/5-floor(n/5))*5-2*double((n/5-floor(n/5))==0)-(n - fix(n))-double(mod(n,1) == 0)-tk);
                delta=1;
                while (value <0.999*Target_effect(i) || value>1.001*Target_effect(i))
                    if value > 1.001*Target_effect(i)
                        n=n-delta;
                        t_trat=7*floor(n/5)+(n/5-floor(n/5))*5-2*double((n/5-floor(n/5))==0)-(n - fix(n))-double(mod(n,1) == 0); 
                        value_integern=floor(n).*alpha.*(2+2.*(a.*2+exp(-a.*2)-1)./(alphabeta.*a^2));
                        d2=2*(n - fix(n));
                        value_fractionaln=alpha.*(d2+2.*(a.*d2+exp(-a.*d2)-1)./(alphabeta.*a^2));
                        value=value_integern+value_fractionaln-lambda*max(0,t_trat-tk);
                        delta=delta/2;
                    elseif value < 0.999*Target_effect(i)
                        n=n+delta;
                        t_trat=7*floor(n/5)+(n/5-floor(n/5))*5-2*double((n/5-floor(n/5))==0)-(n - fix(n))-double(mod(n,1) == 0); 
                        value_integern=floor(n).*alpha.*(2+2.*(a.*2+exp(-a.*2)-1)./(alphabeta.*a^2));
                        d2=2*(n - fix(n));
                        value_fractionaln=alpha.*(d2+2.*(a.*d2+exp(-a.*d2)-1)./(alphabeta.*a^2));
                        value=value_integern+value_fractionaln-lambda*max(0,t_trat-tk);
                        delta=delta/2;
                    end          
                end
                Effect_eq(i)=value;
                n_f(i)=n;
                EQD2(i)=2*n;
            end
        else
    %If there is no issoeffective treatment with 2 Gy/# return NaN
            EQD2=NaN;
            Effect_eq=NaN;
            n_0=NaN;
            n_f=NaN;                
        end  
    end    
end




